function [amplitude, phase_rad] = lockin_amplifier(x, f0, t)
phase_rad = 0;
t_max =fix(max(t)*f0)/f0;
t = t(t<t_max | t==t_max);
x = x(t<t_max | t==t_max);
x = x-mean(x);
for i = 1:5
    % 构造参考信号
    ref_cos = cos(2*pi*f0*t + phase_rad);
    ref_sin = sin(2*pi*f0*t + phase_rad);

    % 正交解调
    I = dot(x, ref_cos)/dot(ref_cos,ref_cos);  % 直流分量 = 实部
    Q = dot(x, ref_sin)/dot(ref_sin,ref_sin);  % 直流分量 = 虚部

    phase_rad = atan2(I,Q)+phase_rad;          % 相位（弧度）
end
amplitude = sqrt(I^2+Q^2);
end
